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ABSTRACT 

The stiffness of stepped bars is determined by calcu- 
lating deflections due to applied loads using the finite 
element method. Nine combinations of step height and fillet 
radius are studied for each of the following cases : 

1. rectangular cross section under axial load; 

2. rectangular cross section under bending load; 

3. circular cross section under axial load; 

4. circular cross section under torsional load. 

A correction parameter, actually a fictitious axial dis- 
placement of the step, is derived. It makes possible ac- 
curate calculation of deflections of stepped bars using 
elementary formulas . 
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I . INTRODUCTION 



There are many mechanical and structural applications 
of bars and beams which have abrupt changes in cross-sec- 
tional area. In the design of these members their deflec- 
tions under expected loads must be considered. If accuracy 
is not essential, the deflections can be determined by 
relatively basic calculations. These calculations involve 
formulas based on the simple stress distributions found in 
constant section members. They do not take into considera- 
tion the reduced stiffness due to stress flow irregularities 
in the vicinity of the step. A stress concentration is said 
to exist at this location. 

Through the numerical analysis technique called the 
finite element method, the deflection of stepped bars or 
beams can be accurately determined. This method is sensi- 
tive to the effect of the stress concentration. 

The purpose of this study was to analyze the overall 
stiffness of stepped bars. An additional goal was to de- 
termine a correction parameter which could be used in de- 
sign to compensate for the reduced stiffness due to the 
stress concenxration . Finixe element analyses were performed 
for stepped bars of various dimensions under axial, bending, 
and torsional loads.- The resulting deflections were deter- 
mined, converted to a common correction parameter, and 
compared. 
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The recent emergence of the finite element method along 
with the digital computer has made feasible the previously 
prohibitive calculations necessary for stiffness analysis. 
The only prior study of stiffness of stepped bars that this 
author could identify was accomplished by F. P. Porter [1]^. 
His study of torsional vibration of irregular shafts in- 
cluded stiffness calculations based on the mathematical 
theory of elasticity using an approximate graphical 
.construction . 



P- 



^Numbers in brackets refer to the List of References, 
80 . 
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II. ELEMENTARY STRESS RELATIONSHIPS FOR THE THREE LOADING 
CASES 



A. BASIC ASSUMPTIONS 

For each of the three loading cases (axial, bending, and 
torsional) , the stepped bar was assumed to be perfectly 
elastic, homogeneous, isotropic, and to exhibit a linear 
stress-strain relationship. In addition, it was assumed 
that the bar was initially straight and all loads were ap- 
plied smoothly to avoid the effect of impact loading. 



B. THE AXIAL LOAD CASE 

A bar subjected to equal and opposite tensile loads will 
deform according to the equation 
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( 1 ) 



where 6 is the elongation in length X, , P is the axial ten- 
sion, A is the cross-sectional area, and E is the modulus 
of elasticity of the bar material. This equation is valid 
provided that the loading is applied uniformly across both 
ends of the bar and that the bar is of constant cross-section. 

For a bar of non-uniform cross section, the elongation 
is estimated by dividing the bar into segments of constant 
cross section and adding the segment elongations. This pro- 
cedure does not take into account the reduced stiffness due 
to the stress concentration located at the change of cross 
section . 
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Figure 1 is a longitudinal section of the bar studied. 



The estimated elongation of the bar is therefore 





] . 



( 2 ) 



and are the respective cross-sectional areas. 

The axial load case was analyzed for both the plane 
stress state and the axisynunetric stress state. In the 
plane stress analysis , the bar was considered to have unit 




Figure 1. Geometry of stepped bar. 

thickness , therefore the cross-sectional area equaled the 
height of the bar. In the axisynunetric stress analysis the 
bar was considered to be of circular cross section. 



c. THE bendi::g case 



For bars subjected to flexure, or pure bending, a mea- 
sure of the beam flexibility is the angle through which the 
ends of the bar rotate. The total angle can be determined 
from the equation 



4 > 



m 

El 



(3) 
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4 



where 4> is the sum of angles <t >2 <t >2 described in Fig. 

2. M is the bending moment and I is the moment of inertia 
about the horizontal centroidal axis of the cross-sectional 




Figure 2. Bending angles. 



area. Equation (3) is derived from the following equation 
which relates the bending moment around the neutral axis of 
an elastic beam to the curvature, <, of the elastic curve. 



K 



M 

El * 



This theory is valid provided: 

1. Plane sections remain plane during bending. 

2. The bar is of constant cross section. 

3. Bending is such that shear strains are negligible. 
As in the axial load case, the stepped bar can be con- 
sidered in two sections, each of constant cross section. 
The estimate for the total bending angle is therefore 




( 4 ) 
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The bar analyzed was of rectangular cross section and 
considered to have unit thickness. The bending moment was 
applied by a linearly varying normal stress across both 
ends of the bar exerting tension on one side of the neutral 
surface and compression on the other side as shown in Figure 
3. 




Figure 3. Bending stress distribution. 

D. THE TORSION CASE 

The equation for angular deflection of a uniform, round, 
cross-section bar subjected to torsion is 



0 




(5) 



9 is the angular deflection, T is the torque, J is the polar 
moment of inertia of the cross-sectional area, and G is the 
modulus of elasticity in shear. It is assumed that plane 
sections perpendicular to the axis of the bar remain plane 
and that a radial line remains straight vihen the bar is 
tv/isted . 

In this study, one end of the bar v/as held fixed and a 
torque v;as applied to the other end such that the shear 
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stress varied linearly from zero- at the centerline to a 
maximum at the outer radius as shown in Fig. 4. 



T 




Figure 4. Torsional stress applied to circular cross section. 

As before, the bar was considered in two sections of con- 
stant cross section and the total angle of twist was the sum 
of the angles for each section. Therefore the estimate for 
the entire bar is 
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III. METHOD OF STIFFNESS COMPARISON 



It is well documente’d u>_riection of a bar of 

non-uniform cross section is slightly greater than would be 
indicated by basic equations (2), (4), and (6). This de- 
creased stiffness is caused by the stress intensification 
at the point of discontinuity. 

One way of accounting for this decreased stiffness is 
to consider the effect of the stress intensification on the 
bar as a whole to be the same as the effect of moving the 
location of the discontinuity. Increasing the length of 
the smaller portion of the bar and decreasing the length of 
the larger portion by the same amount would achieve the same 
effect on its stiffness. 

This apparent change in the location of the bar discon- 
tinuity was calculated and used as. a measure to compare 
stiffness of bars of different dimensions. The size of 
this change was defined as the variable "A" and is illus- 
trated in Fig. 5. 



— H A h — 




Figure 5. Geometric interpretation of stiffness correction 
parameter. 



To calculate A, equations (2), (4), and (6) were altered as 
follows : 




+ ] 



+ - 4 — ] 



5.^ + A 

+ - 4 — ]. 



(7) 

( 8 ) 

(9) 



In each loading case, the deflections required to cal- 
culate 6 , cj) , and 0 were determined by digital computer finite 
element stress analysis. Knowing the load, modulus of elas- 
ticity, and physical dimensions of the bar, the parameter A 
was determined. The A parameter is common to equations (7), 
(8), and (9) and therefore is a means of comparing stiffness 
corrections for the three types of loading. In addition, it 
is independent of the load applied because the magnitude of 
the deflection and the load applied are directly proportional. 

Since the A parameter is a means of correcting for the 
deflecticr. of a s'**Gpped bar under load, it was defined as 
the "stiffness correction parameter. " 
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IV. SOFTWARE 



A. BACKGROUND 

Before the stiffness correction parameter could be de- 
termined, the deflections of the bar under each of the three 
loading conditions had to be determined. Although photo- 
elastic analysis and direct strain measurement are suitable 
for determining the effect of stress concentration, the ideal 
method for overall stiffness evaluation is the finite ele- 
ment stress analysis. 

In very brief and general terms, the finite element 
method is an advanced numerical method that has been made 
possible by the high speed digital computer. When applied 
to structures it can be used to perform two or three di- 
mensional stress analyses. The structure under study is 
divided into subdivisions called finite elements. Adjacent 
elements are connected along the entire inter-element 
boundary. Specific points at regular intervals along the 
element boundaries are designaxea and defined as nodal 
points or nodes. When a load is applied to the structure, 
the displacements of each of these nodal points are calcu- 
lated and used to determine the stresses. 

For each of the problems , maximum use of existing com- 
puter programs was made. For the plane stress analysis of 
■the axial and bending load cases, a program titled PLISOP 
[2] was used. It was developed at the Naval Postgraduate 
School, Monterey, California and has been in use for several 
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years. For the torsional loading case, a program titled 
TORT 2 [3] was used. It was .de.vi?]x>n'=d. at the University of 
Wales, Swansea, and had been iirs Laij-v.;^ locally but had not 
been used or completely debugged. In order to study the 
fourth problem, the axially loaded bar of circular cross 
section, the capability of PLISOP had to be expanded to 
perform an axisyihmetric stress analysis. 

B. MESH GENERATION 

The network of lines describing the division of a struc- 
ture into finite elements is referred to as a mesh. Nor- 
mally a finer mesh more accurately describes the physical 
structure under study and results in a more accurate stress 
analysis. Unfortunately, increasing the number of subdi- 
visions or elements greatly increases the computational ef- 
fort involved in the problem solution. As a result, the 
maximum capabilities of most general purpose digital com- 
puters are rapidly reached as meshes are refined. 

In many stress analyses, the area of interest is only a 
small portion of the structure. In these cases the mesh 
can be very fine in the area of interest and much coarser 
for the rest of the structure. This minimizes the total 
number of elements yet achieves accuracy where it is desired. 
In this study, however, the analysis of the structure as a 
whole was of interest. Therefore it was desirable to obtain 
a mesh as uniformly refined as possible within computer pro- 
gram limitations. The only practical means of maximizing 
accuracy yet minimizing the total number of elements was to 
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take advantage of symmetry. This was done by dividing the 
bar in half along the longitudinal axis. 

It was found satisfactory to take the lengths 5.^, 
the two bar segments equal to the larger height dimension 
D. This provided a nearly undisturbed stress distribution 
at both ends of the bar. Thus the stress intensification 
resulting from the abrupt change of cross section is negli- 
gible at the ends of the bar. 

Based on these considerations, the stepped bar was 
.represented as shown in Fig. 6. In each loading case the 
ratio of d/D and the fillet radius were varied. These vari- 
ations necessitated minor changes in the mesh representation, 
however, this general discussion will refer to the arrange- 
ment of Fig. 6. 




Figure 6. Finite element mesh. 

For this basic geometry, the mesh used consisted of 
twenty-four isoparametric quadrilateral elements. The con- 
cept of the isoparametric element requires lengthy explana- 
tion and is not considered essential to this development. 
The reader is referred to chapter eight of Ref. 4, for the 
background of the isoparametric element. Quadrilateral 
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elements were used as opposed to triangular elements. 

Both of the computer programs used isoparametric quadri- 
lateral elements. Such elements can be linear, quadratic, 
or cubic. In this study, the only curved portion of the 
mesh, the fillet, was adequately described by quadratic 
elements . 

The quadratic, isoparametric, quadrilateral element has 
four corner nodal points, or nodes, and four mid-side nodes. 
The coordinates of each of these nodes are part of the in- 
put data for both the PLISOP and the TORT 2 programs. The 
task of providing this coordinate data as well as data to 
describe how the elements fit together is a considerable 
one. Fortunately, another computer program was developed 
to perform just this task. This program, titled PLIMEG [5], 
generates a mesh as specified. It requires only enough geo- 
metric data to describe the overall area to be discretized. 

The directions for use of PLIMEG are contained in the program. 

The twenty-four element mesh was chosen as an adequate 
representation of the region. A later convergence study, 
v;hich used a finer mesh, confirmed the adequacy of the 
tv;enty-four element mesh for this primarily comparative 
study . 

Illustrations of the meshes used for the various geo- 
metries are shovm in Appendix .A. 

,C. COMPUTER PROGRAM EMPLOYED 

1. The PLISOP program performs a plane stress or a 
plane strain analysis by the finite element method. Detailed 
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directions for its use are contained in the program. The 
bulk of the input data' -hy the mesh generat- 

ing program, PLIMEG , previously ‘mem.i'O'ii’ed . The axial and 
bending effects were achieved by the loads applied and by 
the boundary conditions imposed. 

The load vectors must be derived externally and are ap- 
plied by specifically defining the load and its direction 
at particular nodal points in the mesh. Derivation of the 
force vectors for both axial and bending loading is shown 
in Appendix B, sections 1 and 2, respectively. 

The boundary conditions necessary for the axial loading 
of the bar half-section shown in Fig. 6 consisted of con- 
straining. each node along the bottom surface (bar neutral 
surface) to freedom of movement in the longitudinal direc- 
tion only. In addition, one of the corner nodes along this 
surface was pinned to prevent longitudinal motion as well. 
For the bending load, the nodes along the bottom surface 
were free to move in the vertical direction only. The two 
corner nodes along this surface were pinned to allow only 
rotational displacement at these two points. 

2. The TORT 2 program performs a stress analysis of 
circular sectioned solids with varying diameter using the 
finite element method. The directions for its use are 
contained in a separate booklet accompanying the program. 
The same meshes v;ere used to represent the stepped bar 
except the origin of the coordinate system had to be at 
the lov/er right corner of the bar. The origin for the mesh 
used previously for the PLISOP program v;as at the lower 
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left corner. A relatively simple computer routine can be 
used to renumber the nodal point data provided by the PLIMEG 
program, but the connectivity data was revised by hand. 

Loading is applied by specifically assigning values of 
shear stress at particular nodal points. 

3. To study a bar of circular cross section under axial 
load, an axisymmetric stress analysis capability was re- 
quired. The most direct way of obtaining this capability 
was to alter the PLISOP program. Due to the symmetry of a 
body of revolution, the problem remains mathematically two 
dimensional. The principal difference between the plane 
stress and the axisymmetric stress problems is the intro- 
duction of a fourth component of stress, the "hoop" stress. 
The details of the development of the axisymmetric capability 
are included in Appendix C. 

Since the PLISOP program is essentially unchanged by 
this added capability, the directions for its use are the 
same. Appendix D contains notes on the axisymmetric capa- 
bility intended to amplify the directions contained in the 
program. Appendix E contains a listing of the PLISOP pro- 
gram in its revised form. 

As in the plane stress analyses of the axial and bending 
load cases , the load vectors must also be derived externally 
for the axisymmetric stress analysis. This derivation is 
included in Appendix B, section 3. 
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V. RESULTS 



A. PRESENTATION OF RESULTS 

For each of the loading cases , nine variations of the 
bar dimensions were analyzed. The dimensions and results 
were nondimensionalized by dividing each by the dimension, 

D, which was held constant throughout the study. The 
values of the stiffness correction parameter, A/D, for each 
of the four problems are presented in tabular form. Tables 
I - IV, and graphically. Figs. 7-10. 

From the graphical presentations it is evident that the 
largest correction parameters are required when there is no 
fillet. As the fillet radius is increased, A/D decreases. 

In addition, as the dimension, d, is varied with r held 
constant, the largest correction is required for the ratio 
of d/D equal to one-half. 

The negative values of A/D indicate the bar is stiffer 
than equations (2), (4), and (6) would show. This is ex- 
plained by the fact that these equations, as applied, do 
not take into consideration the increased cross-sectional 
area resulting from the fillet. The fillet area was not 
included in the equations for calculating the correction in 
order to keep the calculations simple and uniform for all 
cases . 

B. COMPARISON OF THE FOUR CASES 

In order to compare the four cases, a plot of A/D versus 
r/D was prepared for each of the three ratios of d/D. Figures 
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Table I 

A/D for Axial Load 
(rectangular cross section) 



\ d/D 








r/DN. 


0 .25 


0.50 


0.75 


0.0 


0 . 19 


0.22 


0.16 


CL. 12.5. 




.0.15. 


0.09 


0.250 


0 . 04 


0.06 


0.0 



Table II 

A/D for Axial Load 
(circular cross section) 



\d/D 

r/D\ 


0.25 


0.50 


0.75 


0 . 0 


0.10 


0.20 


0.12 


0.125 


0.10 


0.13 


0.05 


0 .250 


■D.05 


0.04 


-0.04 



Table III 

A/D for Bending 
(rectangular cross section) 



\d/D 1 


1 


1 1 


1 1 


r/D \ 


1 

0. 25 


: 0.50 1 


0.75 


j 

0.0 i 


o 

o 


! 1 
1 0.12 


1 0 . 11 


0.125 


1 

0 . 0 


1 

0 . 06 


1 ri ^ 4 

0.04 


0.250 


<D 

O 


-0 . 04 


-0.05 



Table IV 

A/D for Torsion 
(circular cross section) 



\d/D 1 


1 


1 


I ! 


r/D\l 


0.25 


0.50 


0 . 75 


0.0 


i 

0.03 


0.06 


0.06 


0.125 


-0.04 


0.0 


1 — J 
o 

o 


0.250 


-0.12 


<D 

O 


^.09 
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Figure 7. Stiffness correction for axial load (rectangular 
cross section). 
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Figure 8. Stiffness correction for axial load (circular 
cross section ) . 
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Figure 9. Stiffness correction for bending (rectangular 
cross section) . 
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Figure 10. Stiffness correction for torsion (circular 
cross section) . 
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11, 12, and 13 show the effect on stiffness of the three 
different types of loading and two types of cross section. 
The most severe type of loading is clearly axial as indi- 
cated by the larger values of A/D. I. M. Allison [6] also 
found the axial loading to be the most severe of the three. 
His conclusions, however, were based on stress concentra- 
tion factors. 

Comparison of the two types of cross section subjected 
to axial load indicates the circular cross section is 
stiffen than the rectangular cross section for a particular 
ratio of d/D. 

C. CONVERGENCE AND UNCERTAINTY 

The normal methods for proving convergence of a finite 
element solution are by further refining the mesh or by 
using higher order elements. Both of these methods were 
investigated. For the basic bar case (D=4, d=2, r=0) under 
axial load, the convergence test results are presented in 
Fig. 14. 

These results show the improvement in accuracy gained 
by using the quadratic rather than the linear elements. Al- 
though it can be assumed that cubic elements would give 
further improvement, other considerations prevailed. Since 
it was desired to cover a wide range of problems , the addi- 
tional computational effort required for cubic elements was 
.deemed an unnecessary luxury. Also, since the geometry of 
the bar v;as relatively simple, consisting of plane or cylin- 
dric surfaces, cubic elements would not provide the same 
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Figure 11. 



Stiffness correction for various loads with 
d/D = .25. 
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Figure 12. Stiffness correction for various loads with 
d/D = .50. 
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Figure 13. 



Stiffness correction for various loads with 
d/D = .75. 
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Figure 14. A/D vs.l/(no. of nodes). 

percentage improvement that the quadratic elements did. 

In order to check the accuracy achieved by mesh refine- 
ment with quadratic elements, the PLISOP program limits were 
altered to accommodate a ninety-six quadratic element mesh. 
Figure 14 indicates this mesh has come very close to con- 
vergence. Based on these results, the uncertainty in the 
values of A/D in Tables I - IV is considered to be the dif- 
ference between the ninety-six element result and the twenty 
four element result, rounded to the second decimal place. 

The uncertainty is therefore +.02. 

An additional test problem v;as run to verify that the 
geometry studied, v;ith dimensions and equal to the 
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larger height dimension, was adequate to ensure that the 
stress intensification at the discontinuity did not signifi- 
cantly affect the stress distribution at the ends of the bar. 
To check this aspect, the length dimensions and ^2 were- 
increased to 1.5 times the larger height dimension. The 
resulting A/D parameter was within .001 of the corresponding 
case with the shorter length, thus confirming the sufficiency 
of the ratios !L^/'D=Si 2 ^^ ~ 1*0. 

The one geometry and loading case which can be compared 
with Porter's work agrees very favorably. For torsion of a 
bar with d/D = 0.75 and r/D = 0.0, the stiffness correction 
parameters agree to the third decimal place. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The stiffness correction parameters obtained are useful 
and applicable in the design of stepped bars. The tables 
and curves included in Section V can be used with confidence 
to the first decimal place of A/D within the range of physi- 
cal bar dimensions studied. 

The results also yield important comparative data for 
stepped bars of various dimensions under the three types of 
loading. The following comparative observations were made: 

1. The value of A/D for stepped bars subjected to axial 
load is larger than the corresponding value for bending or 
torsion . 

2. The stepped bar with ratio' d/D equal to one-half 
yields a larger value of A/D than the other ratios studied. 

3. Fillets reduce stress intensification at the loca- 
tion of the step and therefore increase the stiffness of 
the bar. 

In addition, it was observed that the localized stress 
intensification had a negligible effect on the cross-sectional 
stress distribution at a distance equal to the height of the 
bar. This illustrates the St. Venant principle of rapid 
dissipation of localized stresses. 

B. RECOMMENDATIONS 

It would be of interest to isolate one loading case and 
through use of more refined meshes and cubic elements achieve 
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correction parameters accurate to the third decimal place. 
In addition, a great deal morf^ -coTnpar'r’ti ve information 
could be gathered with a three crxiiicji&icnal finite element 
solution program. Both of these would require altering 
existing computer programs and would greatly increase the 
computer core storage space required for each run in this 
study. 
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APPENDIX A: 



FINITE ELEMENT MESHES EMPLOYED 
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Figure 15. d/D =0.25 r/D 
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Figure 16. d/D = 0.25 r/D = 0.125 
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Figure 17. d/D = 0.25 r/D = 0.250 
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Figure 18. d/D = 0.50 r/D 
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Figure 20. d/D = 0.50 r/D = 0.250 
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Figure 22. d/D =0,75 r/D = 0.125 
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Figure 23. d/D = 0.75 r/D = 0.250 



APPENDIX B: DERIVATION OF EQUIVALENT NODAL FORCES 



1. Plane Stress Analysis of the Axially Loaded Bar 

Boundary stresses and distributed loads are converted 
to equivalent nodal forces by the principle of virtual work. 
For the case of external stress a applied along the boundary 
of an element, the nodal forces are developed as follows: 

6W = [5u^]'^Cf^]® = / a*6u-d(area) 
substituting 

6u = [6u^]'^CN^]^ 

Cf^]® = / aCN^fd(area) (10) 

where 

6W = increment of virtual work 
f. = nodal forces 

X 

a = applied stress 

6u^ = incremental nodal displacements 

6u = the virtual displacement of a typical point 
on the element surface 

N^ = nodal shape functions 

e = superscript, denotes element contribution. 

The quadratic, isoparametric, quadrilateral element 
pictured below was used throughout this study. The shape 
functions are expressed in terms of the local, normalized 
coordinates (C,n). For loading as shovm , the nodal force 
equations become: 
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Figure 24. Quadratic, isoparametric, quadrilateral element. 



dn 



where 

h = height of the element 
d(area) = dy = (h/2)dri for bar of unit thickness. 
Substituting the shape functions and integrating, the result 
for a equal constant is: 
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Where nodes are shared by two adjacent elements, the 
nodal forces are summed. 

2.. Plane Stress Analysis of a Bar Under Bending Load 

The equivalent nodal forces for the linearly varying 
stress v;ere dex'-ived in the same manner as the axially loaded 
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case except for the applied stress. The linearly varying 
stress shown on the element surface, shown below, can be 




Figure 25. Elemental linearly varying stress distribution, 
represented by the equation 






( 12 ) 



Substituting into the equation for nodal forces 




Applying this formula to one end of a particular bar 
spanned by two elements as shown in Fig. 26 and letting ^^=0, 

then = *^3 = and = 0. 

Therefore, the nodal forces along this surface, num- 
bered from top to bottom are: 
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superscript = element 
number 

subscript = element 

node number 



Figure 26. Global bending stress distribution. 
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3. Axisymmetric Stress Analysis of the Axially Loaded Bar 
In this case , the region associated with an element is 
the sector of an annulus with a one radian central angle. 
Figure 27 illustrates the coordinates used with the axisym- 
metric element and the node numbering used in development 
below . 
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Figure 27. Axisymmetric element. 
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Beginning with equation (10), the equivalent nodal forces 
for a uniform external stress are developed as follows: 



[f^]® = /a* [N^]d(area) 
d(area) = r (1 radian) dr = r dr 
r = r^N^ + r 2 N 2 + T 3 N 3 
dr = + r 2 N^ + r 3 Npdn 

where = 3N^/3ri. Then, 

^^3 ■ ’^l 

dr = 2 • 

Therefore , 

d(area) = ^^1^1 ^2^2 ^ ’ 

Substituting back into equation (10), 



iH CM 

4h 4h 

1 


^3'^1 
-a 2 


7(N3_)Mn /N3^N2dn /N3^N3dn" 
/(N2)''dn /N2N3dn 




■’^l' 

^2 


-^3- 




_sym. /(N 3 )^dri_ 




-^3- 



integrating the shape functions 
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APPENDIX C: DEVELOPMENT OF THE AXISYMMETRIC STRESS CAPABILITY 



For an axisymmetrically loaded bar, displacements are 
limited to the axial and radial directions. In this study 
these displacements were identified with the letters u and 
V, respectively. Cylindrical coordinates are applicable 
and are defined in Fig. 28. It can be seen that radial dis- 
placements introduce a circumferential component of strain 



r 




Figure 2 8. Cylindrical coordinates. 



as well’ as a radial component. For problems of this type, 
the general strain-displacement relationships, as developed 
in chapter five of Ref. 4, are 





f 1 


0 U./ ^ 
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^r 




3v/3r 
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y rz _ 




3v/3z + 3u/3r 



The elasticity matrix which links the strains to the stresses 



is derived from the following stress-strain relationships : 




V(0^ + Og)] 

v(o + a„) ] 

2 U 



53 





The resulting elasticity matrix [D] is 




1 v/(l-v) v/(l-v) 0 

1 v/(l-v) 0 



1 



0 

(l-2v)/2(l-v) 



. (17) 



The finite element formulation of the element stiffness 
matrix is carried out as in the plane stress problem, ex- 
cept the integration over the volume of the element must 
take into consideration the fact that the element thickness 
varies with radius. The element stiffness matrix [k] can 
therefore be expressed by the equation, 



where [B] is the matrix derived from the strain-displacement 
relationships and defined such that 



The primary changes to the PLISOP program vjere to in- 
corporate the new [B], [D] , and [k] matrices, as defined 
above. The element shape functions also had to be added, 
since only their derivatives v/ere included PLISOP. 

For further detail on axisymmetric stress analysis by 
the finite element method, the reader is referred to chapter 
five of Ref. 4. 



[k]® = / [B]'^[D][B]r dr dz 



(18) 



[e] = CB][a]- 



0 

where [a] is the element displacement vector. 
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APPENDIX D: NOTES ON AXISYMMETRIC STRESS ANALYSIS CAPABILITY 

OF PLISOP 

1. z and r denote respectively the axial and radial 
coordinates and u and v the corresponding displacements. 

(In PLISOP, z = X and r = y.) 

2. The volume of material associated with an element 
has a thickness equal to the arc subtended by a one radian 
angle vice the body of revolution usually discussed in most 
texts. It turns out the only difference is a factor of 2tt. 

3. Since the axisymmetric analysis involves division 

by r , there cannot be any nodal points with coordinate r=0. 

- 6 

For nodal points along the centerline, substituting r=lxlO 
suffices to prevent the solution from "blowing up," although 
the accuracy of stress values near these nodal points is ad- 
versely affected. 
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APPENDIX E: COMPUTER PROGRAM 
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bars. 



